Self-organization of heterogeneous topology and symmetry breaking in networks with 

adaptive thresholds and rewiring 
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We study an evolutionary algorithm that locally adapts thresholds and wiring in Random Thresh- 
old Networks, based on measurements of a dynamical order parameter. A control parameter p 
determines the probability of threshold adaptations vs. link rewiring. For any p < 1, we find spon- 
taneous symmetry breaking into a new class of self-organized networks, characterized by a much 
higher average connectivity K evo than networks without threshold adaptation (p — 1). While K evo 
and evolved out-degree distributions are independent from p for p < 1, in-degree distributions be- 
come broader when p — * 1, approaching a power- law. In this limit, time scale separation between 
threshold adaptions and rewiring also leads to strong correlations between thresholds and in-degree. 
Finally, evidence is presented that networks converge to self-organized criticality for large N. 

PACS numbers: 05.45.-a, 05.65. +b, 89.75.-k 



Interaction networks in nature often exhibit highly in- 
homogeneous architectures. Examples are scale-free de- 
gree distributions in protein networks [l[ and metabolic 
networks 0] , mostly accompanied by intricate second or- 
der regularities as, for example, community structure 
The emergence of these properties often is explained by 
means of intuitive topology-based models, e.g. prefer- 
ential attachment Q or node duplications Real net- 
works, however, are characterized not only by an evolving 
topology, but also by evolution of function, conveniently 
abstracted in terms of dynamics, i.e. the flow of informa- 
tion or matter on these networks. So far, only few studies 
explicitly consider the more general case of co-evolution 
between network dynamics and -topology 

One example is the question how networks may evolve 
topologies that optimize biologically relevant parameters, 
e.g. flexible adaptation with respect to changing environ- 
ments, or insensitivity against random perturbations of 
topology or dynamics (robustness) (Io| . In this context, 
Kauffman introduced random Boolean networks (RBN) 
to study the dynamics of gene regulatory networks from 
a global perspective 
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It was shown that RBN 
undergo a order-disorder transition at a critical wiring 



density (connectivity) K c = 2 [ll|, |12], [13|, |14|; similar 
results were established for random threshold networks 
(RTN), which constitute a sub-class of RBN [l^, US, H3] ■ 
It has been postulated that evolution should drive dy- 
namical networks towards this 'edge of chaos' to optimize 
adaptive flexibility and robustness [HUH)]. However, no 
mechanism able to generate critically connected networks 
could be provided. 

To address this problem, a RTN-based model was pro- 
posed, linking rewiring of network nodes to local mea- 
surements of a dynamical order parameter, e.g. the aver- 
age activity (magnetization) Q. It was shown that this 
simple, local adaptive mechanism leads to a global self- 
organized critical state in the limit of large system sizes 
N. Subsequently, this principle was generalized to net- 



works of noisy neurons Q and to RBN with evolvable 
logical functions [9(. Interestingly, finite size networks in 
these models evolve a broadly distributed heterogeneous 
in-degree connectivity 0, [l8| . Still, these topological het- 
erogeneities are smaller than those observed in real- world 
networks, presumably because dynamical elements were 
assumed to be homogeneous with respect to their dy- 
namical behavior. While this assumption leads to elegant 
models, it is quite unrealistic, as it becomes apparent e.g. 
in the frequent occurrence of canalizing functions in gene 
regulatory networks, with strong impact on dynamics in 
RBN models [3]. Considering the accumulating experi- 
mental evidence of both close-to criticality [2(| and het- 
erogeneous architecture 2l| in real gene regulatory net- 
works, it is fascinating to speculate about a mechanism 
that might explain both observations: evolution of local 
dynamical heterogeneity and global homeostasis. 

For this purpose, we introduce a minimal model linking 
regulation of activation thresholds and rewiring of net- 
work nodes in RTN to local measurements of a dynamical 
order parameter. A new control parameter p S [0, 1] de- 
termines the probability of rewiring vs. threshold adap- 
tations. We show that the symmetry of the evolutionary 
attractor for p = 1 (no threshold adaptations, rewiring 
only) is broken spontaneously for any p < 1. This 
new universality class of self-organized networks exhibits 
a much higher average connectivity K evo , compared to 
p = 1 networks, however, with a value K evo that is in- 
sensitive to p. In-degree distributions become very broad, 
approaching a flat power-law tail ~ k ir ^^ 4 for p — > 1. Fur- 
ther, we establish the emergence of strong correlations 
between in-degree and thresholds in this limit, while for 
small p, correlations are weak. This indicates that an 
adaptive time-scale separation, with rare events of dy- 
namical diversification and frequent rewiring, can lead to 
emergence of highly inhomogeneous topologies, without 
the need for network growth (as, for example, in prefer- 
ential attachment models). Finally, we present evidence 
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that networks with p < 1 converge to a critical state for 
large N, however, with a finite size scaling significantly 
different from the one found for the case p = 1. 

Dynamics. We consider a network of iV randomly in- 
terconnected binary elements with states o~i = ±1. For 
each site i, its state at time t+1 is a function of the inputs 
it receives from other elements at time t (synchronous 
updates): 



<Ti(t+l) = 



+1 if/i(t)>0 
— 1 else 



with 



JY 



/»(*) = y^CjjQ-jft) + hi 



(1) 



(2) 



The interaction weights cy take discrete values c,j = ±1, 
with Cy = if site i does not receive any input from 
element j. Thresholds hi may vary from node to node, 
taking integer values hi < In the following dis- 

cussion, adaptive changes will be applied to the absolute 
value \hi\, keeping in mind that the sign of hi is always 
negative. 

As a dynamical order parameter, we define the average 
activity A(i) of a site i 



A(i) 



1 



To 



T 2 
1 t=Ti 



<Ji(t). 



(3) 



Notice that a frozen site, i.e. a site that does not change 
its state, has |A(z)| = 1, whereas an active site has 
\A(i)\ < 1. 




Ai< 1 (active) A, = 7 (frozen) 



FIG. 1: ie/t: with probability p, active nodes loose one of 
their inputs, with probability 1 — p they increase their (ab- 
solute) threshold \hi\. Right: frozen nodes show the opposite 
behavior. 

Topology evolution. Let us now discuss a particular 
evolutionary scheme that couples local adaptations of 
both the number of inputs and of thresholds to a site's 
average activity. Analyzing Eq. (1) and Eq. (J2]), we 
realize that the activity of a site i can be controlled in 
two ways: if i is frozen, it can increase the probability to 
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FIG. 2: Upper panel: Evolution of the average connectivity 
K of threshold networks, using the adaptive algorithm (cf. 
Fig. 1), for N = 512 and initial connectivity Kini = 1. Time 
series for five different values of p are shown. Lower panel: 
The same for the average threshold h. 
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FIG. 3: Upper four curves: Evolutionary mean values K evo of 
the average connectivity, as a function of p; system sizes from 
top to bottom: N = 512, N = 256, N = 128 and N = 64. 
Lower four curves: The same for the evolutionary mean values 
|/iet)o|of the average absolute threshold. Statistics was taken 
over 10 6 evolutionary steps, after a transient of 4 • 10 6 steps. 



change its state by either increasing its number of inputs 
ki — ► ki + 1, or by making its threshold hi < less neg- 
ative, i.e. \hi\ — y \hi \ — 1. If i is active, it can reduce its 
activity by adapting either ki — > ki — 1 or \hi\ — » \hi \ + 1. 
This adaptive scheme is realized in the following algo- 
rithm (see also Fig. 1): 

1. Create a random network with average connectivity 



Ki n i > and average threshold hi 



0. 
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FIG. 4: Line-pointed curves: in-degree distributions of 
evolved networks, data points only: the corresponding out- 
degree distributions ((A) p = 0.3, (+) p = 0.5, (x) p = 0.8, 
(*) p = 0.95, (□) p = 0.99). Statistics was gathered over 10 6 
evolutionary steps, after a transient of 4 • 10 6 steps. Networks 
had size N = 512. The dashed line has slope —3/4. 




FIG. 5: Average number {ki„) of inputs for a given node in 
evolving networks, as a function of the respective nodes (ab- 
solute) threshold \h\. Statistics was taken over 10 6 rewiring 
steps, after a transient of 4 ■ 10 6 steps. For all values p < 1, a 
clear positive correlation between ki n and \h\ is found. 



2. Select a random initial state Q% n % = (cri, ...,<7jv). 

3. Iterate network dynamics for T timesteps. 

4. Select a network site i at random and measure its av- 
erage activity At over the last T/2 updates. 

5. Adapt fcj and hi in the following way: 

- If \Ai\ < 1, then ki—>ki — l with probability p (removal 
of one randomly selected input). With probability 1 — p, 
adapt \ht\ — ► \h% \ + 1 instead. 

- If | A; | = 1, then fcj — > ki + 1 with probability p (addi- 
tion of a new input from a randomly selected site). With 
probability 1 —p, adapt \h%\ — > — 1 instead. If hi = 0, 
let its value unchanged. 

6. Go back to step 3. 

If the control parameter p takes values p > 1/2, 
rewiring of nodes is favored, whereas for p < 1/2 thresh- 
old adaptations are more likely. Notice that the model 
introduced in 0] is contained as the limiting case p = 1 



FIG. 6: Average fraction y(N) of damaged nodes, 200 updates 
after a one-bit perturbation at time t = 0, for different p, as 
a function of system size N. The lined curve is a fit of the 
average scaling behavior. 



(rewiring only and hi = const. = for all sites). 

Results. After a large number of adaptive cycles, 
networks self-organize into a global evolutionary steady 
state. An example is shown in Figure 2 for networks 
with N = 512: starting from an initial value Ki n i = 1, 
the networks' average connectivity K first increases, and 
then saturates around a stationary mean value K evo ; 
similar observations are made for the average threshold 
h (Fig. 2, lower panel). The non-equilibrium nature 
of the system manifests itself in limited fluctuations of 
both K and h around K evo and h evo . Regarding the de- 
pendence of K with respect to p, we make the interest- 
ing observation that it changes non-monotonically. Two 
cases can be distinguished: when p = 1, K stabilizes 
at a very sparse mean value K evo , e.g. for N — 512 at 
K evo = 2.664 ± 0.005. When p < 1, the symmetry of 
this evolutionary steady state is broken. Now, K con- 
verges to a much higher mean value K evo w 43.5 ± 0.3 
(for N = 512), however, the particular value which is fi- 
nally reached is independent of p. The latter observation 
is made rigorous from measurements of K evo for differ- 
ent N over 10 6 evolutionary steps, after systems have 
reached the steady state. While K evo obviously depends 
on the system size N, curves are very flat with respect to 
p (Fig. 3, upper four curves); the same holds for \h evo \ 
(Fig. 3, lower four curves). On the other hand, conver- 
gence times T con needed to reach the steady state are 
strongly influenced by p: T con (p) diverges when p ap- 
proaches 1 (compare Fig. 2 for p = 0.99). We conclude 
that p determines the adaptive time scale. This is also 
reflected by the stationary in-degree distributions p(ki n ) 
that vary considerably with p (Fig. 4); when p — > 1, 
these distributions become very broad. The numerical 
data suggest that a power law 



limp(fc m ) oc /c in 7 



(4) 



with 7 « 3/4 ± 0.03 is approached in this limit (cf. Fig. 
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4, dashed line). At the same time, it is interesting to 
notice that the evolved out-degree distributions are much 
narrower and completely insensitive to p (Fig. 4, data 
points without lines). 

How can one understand the emergence of broad in- 
degree distributions for with increasing pi Evidently, life 
times of both low thresholds \hi\ and high thresholds 
\hi\ ^S> become long for p — ► 1. Since sites with low 
thresholds tend to be active and hence, on average, loose 
links, while sites with high thresholds tend to freeze and 
hence, on average, aquire new links, we would indeed 
expect that p(ki n ) is broadened for p — > 1. On the other 
hand, for p — > 0, frequent adaptive changes of thresholds 
prevent long sequences of both frozen or highly active 
states, and hence make emergence of strong local wiring 
heterogeneities less probable. If this idea is correct, we 
would expect that, in the limit p — > 1, the in-degree of 
sites should exhibit a strong positive correlation to their 
thresholds, while for p — > these correlations should be 
less pronounced. This is indeed exactly what we observe. 
For p — 0.99, the average in-degree (k in ) of a given node, 
as a function of its threshold \h\, shows a steep increase, 
while the corresponding curve is relatively flat for p = 0.3 
(Fig. 5). 

An interesting question is whether the networks with 
p < 1 still approach a self-organized critical state for large 
N, as it was found for the case p = 1 0- Since networks 
now evolve more densely wired, non-trivial topologies, 
this question has to be answered by application of a dy- 
namical criterion. For this purpose, we studied damage 
spreading: after each adaptive step, dynamics was run 
from an initial system state a and again from a direct 
neighbor state a 1 differing in one bit; after t = 200 up- 
dates, the Hamming distance d between both trajectories 
was measured and the average fraction of damaged nodes 
y(t) =: d/N was determined. Figure 6 shows y, averaged 
over 10 6 evolutionary steps, as a function of N . We find 
that the finite networks investigated here are all super- 
critical, however, y decreases monotonically with N . The 
average scaling behavior can be fit by 



y(N) 



[ln(iV)]- 



(5) 



with o = 0.77 ± 0.02 and (3 = 0.917 ± 0.01. This de- 
pendence indicates that y — 0, i.e. the critical transi- 
tion form chaotic to frozen dynamics, is approached for 
large N . Notice, however, that convergence is logarith- 
mic, whereas for p — 1 power laws were found [3, 
Again, this indicates that p < 1 networks constitute an 
entirely new universality class. 

To summarize, we studied a model of network evolu- 
tion that couples both rewiring of inputs and adaptation 
of activation thresholds to local measurements of a dy- 
namical order parameter. A control parameter p deter- 
mines the probability of threshold adaptations vs. link 
rewiring. While for p — 1 (rewiring only, no thresh- 
old adapttation) networks evolve a self-organized criti- 



cal state with a sparse average connectivity K evo m 2, 
for any p < 1 (both rewiring and threshold adaptation) 
networks evolve a significantly more dense wiring, with 
broad heterogeneous in-degree distributions approaching 
a power-law ~ k ir ^^ 4 for p — ► 1. In this limit, time scale 
separation between rare threshold adaptations and fre- 
quent rewiring leads to emergence of strong correlations 
between thresholds and in-degree. We presented evidence 
that, in the limit of large N, networks logarithmically ap- 
proach a self-organized critical state. 

Our model presents a novel mechanism leading to 
co-evolution of topological and dynamical heterogeneity 
with robust homeostatic regulation, the latter reflected 
e.g. by the insensitivity of the evolved average connec- 
tivity with respect to p. Since similar - seemingly contra- 
dicting - observations are also made in experimental data 
of, e.g., gene regulatory networks 2(J 21 1, it is interesting 
to speculate that similar mechanisms might be at work 
in the evolution of biological networks. 
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